ERASE SOURCE.READ$$
FTN READ$$,N,ISD=N
      SUBROUTINE READ (I, A, NA, B, NB, C, NCX, D, ND, E, NE)
      DIMENSION A(1), B(1), C(1), D(1), E(1)
      DIMENSION NA(2), NB(2),NCX(2), ND(2), NE(2),NZ(2)
      DOUBLE PRECISION A, B, C, D, E
      READ(5,100) LAB,            NZ(1), NZ(2)
      CALL READ1(A, NA,NZ,  LAB)
      IF(I .EQ. 1) GO TO 999
      READ(5,100) LAB,            NZ(1), NZ(2)
      CALL READ1(B, NB,NZ,  LAB)
      IF(I .EQ. 2) GO TO 999
      READ(5,100) LAB,            NZ(1), NZ(2)
      CALL READ1(C, NCX,NZ, LAB)
      IF(I .EQ. 3) GO TO 999
      READ(5,100) LAB,            NZ(1), NZ(2)
      CALL READ1(D, ND,NZ,  LAB)
      IF(I .EQ. 4) GO TO 999
      READ(5,100) LAB,            NZ(1), NZ(2)
      CALL READ1(E, NE,NZ,  LAB)
  100 FORMAT(A4,4X,2I4)
  999 RETURN
      END
ERASE SOURCE.RDTITL$$
FTN RDTITL$$,N,ISD=N
      SUBROUTINE RDTITL
      COMMON /LINES/NLP,LIN,TITLE(23)
      READ  (5,100)  (TITLE(I),I=1,18)
  100 FORMAT  (18A4)
      CALL LNCNT(100)
      RETURN
      END
ERASE SOURCE.PRNT$$
FTN PRNT$$,N,ISD=N
      SUBROUTINE  PRNT(AR,NAR,NAM,IP)
C   SUBR PRNT   PRINTS DOUBLE PRECISION MATRIX
      COMMON  /FORM/NEPR,FMT1(6),FMT2(6)
      COMMON/LINES/NLP,LIN,TITLE(23)
      COMMON  /MAX/MAXRC
C- NOTE NLP NO. LINES/PAGE VARIES WITH THE INSTALLATION.
      DATA KZ,KW,KB /1H0,1H1,1H /
      REAL*8 AR                                                         D
      DIMENSION AR(1),NAR(2)
      NAME = NAM
C-IF IP =1,HEADLINE SAME PAGE, IF IP =2,  HEADLINE, NEW PAGE
C     IP=3,  NO HEADLINE, SAME PAGE,    IP=4, NO HEADLINE, NEW PAGE
      II = IP
      NR=NAR(1)
      NC=NAR(2)
      NLST = NR * NC
      IF(NLST .GT. MAXRC .OR. NLST .LT. 1.OR.NR.LT.1)  GO TO 16
      IF(NAME  .EQ. 0) NAME = KB
C- SKIP HEADLINE IF REQUESTED.
      GO TO (11,10,132,12),       II
   10 CALL LNCNT(100)
   11 CALL LNCNT(2)
    3 WRITE(6,177) KZ,NAME,NR,NC
177   FORMAT(A1,5X,A4,8H  MATRIX,5X,I3,5H ROWS,5X,I3,8H COLUMNS)
      GO TO 13
   12 CALL LNCNT(100)
      GO TO 13
  132 CALL LNCNT(2)
      WRITE (6,891)
  891 FORMAT (1H0)
C- BELOW COMPUTE NR OF LINES/ ROW --DECIDE IF 1 EXTRA BLANK LINE
   13 J=(NC-1)/NEPR+1
C- WHY ALWAYS ADD 1 LINE- BECAUSE IF MULTIPLE, USE 1 BLK LINE EXTRA.
      NLPW = J
      JST = 1
C- COMPUTE LAST ROW POSITION -1 BELOW
      NLST = NLST -NR
      MN=NC
      IF  (NC.GT.NEPR)   MN=NEPR
      KLST=NR*(MN-1)
91    CONTINUE
      DO 912 J = JST, NR
      CALL LNCNT(NLPW)
      KLST = KLST +1
      WRITE (6,FMT1)    (AR(N), N = J, KLST, NR)
      IF  (NC.LE.NEPR)   GO TO 912
      NLST = NLST +1
      KNR=KLST+NR
      WRITE (6,FMT2)(AR(N),N=KNR,NLST,NR)
912   CONTINUE
      RETURN
   16 CALL LNCNT(1)
      WRITE  (6,916)   NAM,NAR
  916 FORMAT  (' ERROR IN PRNT  MATRIX 'A4,' HAS NA=',2I6)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.LNCNT$$
FTN LNCNT$$,N,ISD=N
      SUBROUTINE LNCNT  (N)
      COMMON/LINES/NLP,LIN,TITLE(23)
      LIN=LIN+N
      IF  (LIN.LE.NLP)   GO TO 20
      WRITE  (6,1010)  (TITLE(I),I=1,23)
 1010 FORMAT  (1H1,23A4/)
      LIN=2+N
      IF  (N.GT.NLP)  LIN=2
   20 RETURN
      END
ERASE SOURCE.ADD$$
FTN ADD$$,N,ISD=N
      SUBROUTINE ADD (A,NA,B,NB,C,NC)
      DIMENSION A(1),B(1),C(1),NA(2),NB(2),NC(2)
      COMMON  /MAX/MAXRC
      DOUBLE PRECISION A,B,C
      IF((NA(1).NE.NB(1)).OR.(NA(2).NE.NB(2))) GO TO 999
      NC(1)=NA(1)
      NC(2)=NA(2)
      L=NA(1)*NA(2)
      IF  (NA(1).LT.1.OR.L.LT.1.OR.L.GT.MAXRC)  GO TO 999
      DO 300 I=1,L
  300 C(I)=A(I)+B(I)
      GO TO 1000
  999 CALL LNCNT (1)
      WRITE(6,50) NA,NB
   50 FORMAT  (' DIMENSION ERROR IN ADD     NA='2I6,5X,'NB='2I6)
      CALL ASPERR
 1000 RETURN
      END
ERASE SOURCE.SUBT$$
FTN SUBT$$,N,ISD=N
      SUBROUTINE SUBT(A,NA,B,NB,C,NC)
      DIMENSION A(1),B(1),C(1),NA(2),NB(2),NC(2)
      DOUBLE PRECISION A,B,C
      COMMON  /MAX/MAXRC
      IF((NA(1).NE.NB(1)).OR.(NA(2).NE.NB(2))) GO TO 999
      NC(1)=NA(1)
      NC(2)=NA(2)
      L=NA(1)*NA(2)
      IF  (NA(1).LT.1.OR.L.LT.1.OR.L.GT.MAXRC)  GO TO 999
      DO 300 I=1,L
  300 C(I)=A(I)-B(I)
      GO TO 1000
  999 CALL LNCNT (1)
      WRITE(6,50) NA,NB
   50 FORMAT  (' DIMENSION ERROR IN SUBT    NA='2I6,5X,'NB='2I6)
      CALL ASPERR
 1000 RETURN
      END
ERASE SOURCE.MULT$$
FTN MULT$$,N,ISD=N
      SUBROUTINE MULT(A,NA,B,NB,C,NC)
      DIMENSION A(1),B(1),C(1),NA(2),NB(2),NC(2)
      DOUBLE PRECISION A,B,C
      COMMON  /MAX/MAXRC
      NC(1) = NA(1)
      NC(2) = NB(2)
      IF(NA(2).NE.NB(1)) GO TO 999
      NAR = NA(1)
      NAC = NA(2)
      NBC = NB(2)
      NAA=NAR*NAC
      NBB=NAR*NBC
      IF  (NAR.LT.1.OR.NAA.LT.1.OR.NAA.GT.MAXRC.OR.NBB.LT.1.OR.
     1    NBB.GT.MAXRC)  GO TO 999
      IR = 0
      IK=-NAC
      DO 300 K=1,NBC
      IK = IK + NAC
      DO 300 J=1,NAR
      IR=IR+1
      IB=IK
      JI=J-NAR
      C(IR)=0.
      DO 300 I=1,NAC
      JI = JI + NAR
      IB=IB+1
      C(IR)=C(IR)+A(JI)*B(IB)
  300 CONTINUE
      GO TO 1000
  999 CALL LNCNT (1)
      WRITE(6,500) (NA(I),I=1,2),(NB(I),I=1,2)
  500 FORMAT  (' DIMENSION ERROR IN MULT    NA='2I6,5X,'NB='2I6)
      CALL ASPERR
 1000 RETURN
      END
ERASE SOURCE.SCALE$$
FTN SCALE$$,N,ISD=N
      SUBROUTINE SCALE (A, NA, B, NB, S)
      DIMENSION A(1),B(1),NA(2),NB(2)
      COMMON  /MAX/MAXRC
      DOUBLE PRECISION A, B, S
      NB(1) = NA(1)
      NB(2) =NA(2)
      L = NA(1)*NA(2)
      IF  (NA(1).LT.1.OR.L.LT.1.OR.L.GT.MAXRC)  GO TO 999
      DO 300 I=1,L
  300 B(I)=A(I)*S
 1000 RETURN
  999 CALL LNCNT(1)
      WRITE  (6,50) NA
   50 FORMAT  (' DIMENSION ERROR IN SCALE   NA='2I6)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.TRANP$$
FTN TRANP$$,N,ISD=N
      SUBROUTINE TRANP(A,NA,B,NB)
      DIMENSION A(1),B(1),NA(2),NB(2)
      DOUBLE PRECISION A,B
      COMMON  /MAX/MAXRC
      NB(1)=NA(2)
      NB(2)=NA(1)
      NR=NA(1)
      NC=NA(2)
      L=NR*NC
      IF  (NR   .LT.1.OR.L.LT.1.OR.L.GT.MAXRC)  GO TO 999
      IR=0
      DO 300 I=1,NR
      IJ=I-NR
      DO 300 J=1,NC
      IJ=IJ+NR
      IR=IR+1
  300 B(IR)=A(IJ)
      RETURN
  999 CALL LNCNT(1)
      WRITE  (6,50)  NA
   50 FORMAT  (' DIMENSION ERROR IN TRANP   NA='2I6)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.INV$$
FTN INV$$,N,ISD=N
      SUBROUTINE INV(A,NA,DET,L)
      DIMENSION A(1), L(1),NA(2)
      DOUBLE PRECISION A, DET, BIGA ,  HOLD
      COMMON  /MAX/MAXRC
      IF  (NA(1).NE.NA(2))  GO TO 600
C     SEARCH FOR LARGEST ELEMENT
      DET= 1.
      N=NA(1)
      NSQ=N*N
      IF  (N.LT.1.OR.NSQ.GT.MAXRC)   GO TO 600
      NK = - N
      DO 80 K= 1, N
      NK = NK + N
      L(K) = K
      NPK=N+K
      L(NPK)=K
      KK = NK + K
      BIGA = A(KK)
      DO 20 J= K, N
      IZ = N*(J - 1)
      DO   20 I= K, N
      IJ = IZ + I
   10 IF(DABS(BIGA) - DABS(A(IJ))) 15, 20, 20
   15 BIGA = A(IJ)
      L(K) = I
      NPK=N+K
      L(NPK)=J
   20 CONTINUE
C     INTERCHANGE ROWS
      J = L(K)
      IF(J - K) 35, 35, 25
   25 KI = K - N
      DO 30 I = 1, N
      KI = KI + N
      HOLD = -A(KI)
      JI = KI - K + J
      A(KI) = A(JI)
   30 A(JI) = HOLD
C     INTERCHANGE COLUMNS
   35 NPK=N+K
      I=L(NPK)
      IF (I - K) 45, 45, 38
   38 JP = N*(I - 1)
      DO 40 J= 1, N
      JK = NK + J
      JI = JP + J
      HOLD = -A(JK)
      A(JK) = A(JI)
   40 A(JI) = HOLD
C     DIVIDE COLUMN BY MINUS PIVOT(VALUE OF PIVOT ELEMENTS IS CONTAINED
C     IN BIGA)
   45 IF(BIGA) 48, 46, 48
   46 DET = 0.
      CALL LNCNT (1)
      KKK=KK-1
      WRITE  (6,1046)  KKK
 1046 FORMAT  (' INV ERROR  DETERMINANT OF A=0   RANK OF A=',I4)
      CALL ASPERR
      RETURN
   48 DO 55 I= 1, N
      IF (I - K) 50, 55, 50
   50 IK = NK + I
      A(IK) =-A(IK)/(BIGA)
   55 CONTINUE
C     REDUCE MATRIX
      DO 65 I= 1, N
      IK = NK + I
      HOLD = A(IK)
      IJ = I - N
      DO 65 J= 1, N
      IJ = IJ + N
      IF(I  - K ) 60, 65, 60
   60 IF(J- K) 62, 65, 62
   62 KJ = IJ -I + K
      A(IJ) = HOLD* A(KJ) + A(IJ)
   65 CONTINUE
C     DIVIDE ROW BY PIVOT
      KJ = K - N
      DO 75 J= 1, N
      KJ = KJ + N
      IF(J - K) 70, 75, 70
   70 A(KJ) = A(KJ)/BIGA
   75 CONTINUE
C     PRODUCT OF PIVOTS
      DET=DET*BIGA
C     REPLACE PIVOT BY RECIPROCAL
      A(KK) = 1./BIGA
   80 CONTINUE
C     FINAL ROW AND COLUMN INTERCHANGE
      K = N
  100 K = K - 1
      IF(K) 150, 150, 105
  105 I = L(K)
      IF (I - K ) 120, 120, 108
  108 JQ = N*( K - 1)
      JR = N*(I- 1)
      DO 110 J= 1, N
      JK = JQ + J
      HOLD = A(JK)
      JI = JR + J
      A(JK) = - A(JI)
  110 A(JI) = HOLD
  120 NPK=N+K
      J=L(NPK)
      IF(J - K) 100, 100, 125
  125 KI = K - N
      DO 130 I= 1, N
      KI = KI + N
      HOLD = A(KI)
      JI = KI - K + J
      A(KI) = - A(JI)
  130 A(JI) = HOLD
      GO TO 100
  150 RETURN
  600 CALL LNCNT (1)
      WRITE  (6,1600)  NA
 1600 FORMAT  (' INV REQUIRES SQUARE MATRIX   NA=',2I4)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.NORM$$
FTN NORM$$,N,ISD=N
      SUBROUTINE NORM(A,NA,ANORM)
      DIMENSION NA(2),A(1)
      DOUBLE PRECISION A,ANORM,SUM,ROWMAX,COLMAX
      COMMON  /MAX/MAXRC
      NAR = NA(1)
      NAC = NA(2)
      L=NAR*NAC
      IF  (NAR  .LT.1.OR.L.LT.1.OR.L.GT.MAXRC)  GO TO 999
      COLMAX = 0.
      ROWMAX = 0.
      K = 0
      DO 300 I = 1,NAC
      SUM = 0.
      DO 301 J = 1,NAR
      K = K + 1
  301 SUM = SUM + DABS(A(K))
      IF (COLMAX.LT.SUM) COLMAX = SUM
  300 CONTINUE
      DO 302 I = 1,NAR
      SUM = 0.
      K = I - NAR
      DO 303 J = 1,NAC
      K = K + NAR
  303 SUM = SUM + DABS(A(K))
      IF (ROWMAX.LT.SUM) ROWMAX=SUM
  302 CONTINUE
      ANORM = DMIN1(COLMAX,ROWMAX)
      RETURN
  999 CALL LNCNT (1)
      WRITE  (6,50)  NA
   50 FORMAT  (' DIMENSION ERROR IN NORM    NA='2I6)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.UNITY$$
FTN UNITY$$,N,ISD=N
      SUBROUTINE UNITY(A,NA)
      DIMENSION A(1),NA(2)
      DOUBLE PRECISION A
      IF(NA(1).NE.NA(2)) GO TO 999
      CALL SCALE(A,NA,A,NA,0.D0)
      J = - NA(1)
      NAX = NA(1)
      DO 300 I=1,NAX
      J=NAX +J+1
  300 A(J)=1.
      GO TO 1000
  999 CALL LNCNT (1)
      WRITE(6, 50)(NA(I),I=1,2)
   50 FORMAT  (' DIMENSION ERROR IN UNITY   NA='2I6)
      CALL ASPERR
 1000 RETURN
      END
ERASE SOURCE.TRCE$$
FTN TRCE$$,N,ISD=N
      SUBROUTINE TRCE  (A,NA,TR)
      DOUBLE PRECISION A(1),TR
      DIMENSION NA(2)
      COMMON  /MAX/MAXRC
      IF (NA(1).NE.NA(2)) GO TO 600
      TR=0.D0
      N=NA(1)
      NN=N*N
      IF  (N.LT.1.OR.NN.GT.MAXRC)   GO TO 600
      DO 10 I=1,N
      M=I+N*(I-1)
   10 TR=TR+A(M)
      RETURN
  600 CALL LNCNT(1)
      WRITE (6,1600) NA
 1600 FORMAT (' TRACE REQUIRES SQUARE MATRIX    NA=',2I6)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.EQUATE$$
FTN EQUATE$$,N,ISD=N
      SUBROUTINE EQUATE(A,NA,B,NB)
      DIMENSION A(1),B(1),NA(2),NB(2)
      DOUBLE PRECISION A, B
      COMMON  /MAX/MAXRC
      NB(1) = NA(1)
      NB(2) =NA(2)
      L=NA(1)*NA(2)
      IF  (NA(1).LT.1.OR.L.LT.1.OR.L.GT.MAXRC)  GO TO 999
      DO 300 I=1,L
  300 B(I)=A(I)
 1000 RETURN
  999 CALL LNCNT (1)
      WRITE  (6,50)  NA
   50 FORMAT  (' DIMENSION ERROR IN EQUATE  NA='2I6)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.JUXTC$$
FTN JUXTC$$,N,ISD=N
      SUBROUTINE JUXTC(A,NA,B,NB,C,NC)
      DIMENSION A(1),B(1),C(1),NA(2),NB(2),NC(2)
      DOUBLE PRECISION A,B,C
      COMMON  /MAX/MAXRC
      IF (NA(1).NE.NB(1)) GO TO 600
      NC(1)=NA(1)
      NC(2)=NA(2)+NB(2)
      L=NA(1)*NA(2)
      NNC=NC(1)*NC(2)
      IF  (NA(1).LT.1.OR.L.LT.1.OR.L.GT.MAXRC)  GO TO 600
      IF  (NC(2).LT.1.OR.NNC.GT.MAXRC)   GO TO 600
      MS=NA(1)*NA(2)
      DO 10 I=1,MS
   10 C(I)=A(I)
      MBS=NA(1)*NB(2)
      DO 20 I=1,MBS
      J=MS+I
   20 C(J)=B(I)
      RETURN
  600 CALL LNCNT(1)
      WRITE (6,1600) NA,NB
 1600 FORMAT (' DIMENSION ERROR IN JUXTC,  NA=',2I6,5X,'NB=',2I6)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.JUXTR$$
FTN JUXTR$$,N,ISD=N
      SUBROUTINE JUXTR(A,NA,B,NB,C,NC)
      DIMENSION A(1),B(1),C(1),NA(2),NB(2),NC(2)
      DOUBLE PRECISION A,B,C
      COMMON  /MAX/MAXRC
      IF(NA(2).NE.NB(2))GO TO 600
      NC(2)=NA(2)
      NC(1)=NA(1)+NB(1)
      L=NA(1)*NA(2)
      NNC=NC(1)*NC(2)
      IF  (NA(1).LT.1.OR.L.LT.1.OR.L.GT.MAXRC)  GO TO 600
      IF  (NC(2).LT.1.OR.NNC.GT.MAXRC)   GO TO 600
      MCA=NA(2)
      MRA=NA(1)
      MRB=NB(1)
      MRC=NC(1)
      DO 10 I=1,MCA
      DO 10 J=1,MRA
      K=J+MRA*(I-1)
      L=J+MRC*(I-1)
   10 C(L)=A(K)
      DO 20 I=1,MCA
      DO 20 J=1,MRB
      K=J+MRB*(I-1)
      L=MRA+J+MRC*(I-1)
   20 C(L)=B(K)
      RETURN
  600 CALL LNCNT(1)
      WRITE(6,1600) NA,NB
 1600 FORMAT(' DIMENSION ERROR IN JUXTR,  NA=',2I6,5X,'NB=',2I6)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.EAT$$
FTN EAT$$,N,ISD=N
      SUBROUTINE EAT(A, NA, TT, B, NB, C, NC, DUMMY, KDUM)
      DIMENSION A(1),B(1),DUMMY(1),NA(2),NB(2),ND(2),C(1),NC(2)
      DOUBLE PRECISION A, T, TT, ANAA, TMAX, B, DUMMY, C , S, SC
      COMMON  /MAX/MAXRC
      NR=NA(1)
      NCC=NA(2)
      NC(1)=NR
      NC(2)=NCC
      NB(1)=NR
      NB(2)=NCC
      LD=NR*NCC
      IF (NR.NE.NCC.OR.NR.LT.1           .OR.LD.GT.MAXRC) GO TO 998
      NDD=2*NA(1)*NA(1)
      IF(KDUM  .LT.NDD) GO TO 998
      NDD= NA(1)*NA(1)+1
      T=TT
      CALL NORM(A,NA,ANAA)
      TMAX= 10.01/ANAA
      K=0
  101 IF (TMAX-T ) 103,104,104
  103 K=K+1
      T=TT/2**K
      IF (K-1000) 101,102,102
  104 SC=T
      CALL SCALE(A,NA,A,NA,T)
      DO 401 I= 1, 2
  401 NB(I) = NA(I)
      CALL UNITY(B,NB)
      CALL SCALE(B,NB,DUMMY(1),NB,T)
      S = T/2.
      CALL SCALE(A,NA,DUMMY(NDD),NA,S)
      N = 35
      II=2
      CALL ADD (DUMMY(1),NA,DUMMY(NDD),NA,DUMMY(NDD),NA)
      CALL ADD(A,NA,B,NB,DUMMY(1),NA)
      CALL EQUATE(A,NA,C,NC)
  106 CALL MULT(A,NA,C,NC,B,NB)
      S=1.D0/II
      CALL SCALE(B,NB,C,NC,S)
      S = T/(II+1)
      CALL SCALE(C,NC,B,NB,S)
      CALL ADD(B,NB,DUMMY(NDD),NB,DUMMY(NDD),NB)
      CALL ADD(C,NC,DUMMY(1),NC,DUMMY,NC)
      N=N-1
      IF (N) 107,107,105
  105 II=II+1
      GO TO 106
  107 CALL EQUATE(DUMMY(1), NB, B, NB)
      IF (K) 109,108,212
  109 CALL LNCNT (1)
      WRITE (6,110)
  110 FORMAT (' ERROR IN EAT    K IS NEGATIVE')
  112 IF (K-1) 213,212,212
  213 K=1
  212 DO 111 J=1,K
      T=2*T
      CALL EQUATE(B, NB, DUMMY(1), NB)
      CALL MULT(DUMMY(1), NA, DUMMY(NDD), NA,C,NC)
      CALL ADD(DUMMY(NDD), NC, C, NC, DUMMY(NDD), NC )
  111 CALL MULT(DUMMY(1), NB, DUMMY(1), NB, B, NB)
      TT=T
  108 CONTINUE
      CALL EQUATE(DUMMY(NDD),NC,C,NC)
      S=1.D0/SC
      CALL SCALE(A,NA,A,NA,S)
      RETURN
  102 CALL LNCNT (1)
      WRITE (6,119)
  119 FORMAT  (' ERROR IN EAT    K=1000')
      CALL ASPERR
      RETURN
  998 CALL LNCNT (1)
      WRITE(6, 50)  NA,KDUM,NDD
   50 FORMAT  (' DIMENSION ERROR IN EAT     NA='2I6,   'KDUM='I5,5X,
     1   'KDUM(MIN)='I5)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.ETPHI$$
FTN ETPHI$$,N,ISD=N
      SUBROUTINE ETPHI(A,NA,TT,B,NB,DUMMY,KDUM)
      DIMENSION A(1),B(1),DUMMY(1),NA(2),NB(2),ND(2)
      DOUBLE PRECISION A, T, TT, ANAA, TMAX, B, DUMMY ,S, SC
      COMMON  /MAX/MAXRC
      NR=NA(1)
      NCC=NA(2)
      NB(1)=NR
      NB(2)=NCC
      LD=NR*NCC
      IF (NR.NE.NCC.OR.NR.LT.1           .OR.LD.GT.MAXRC) GO TO 998
      NDD=2*NA(1)*NA(1)
      IF(KDUM  .LT.NDD) GO TO 998
      NDD= NA(1)*NA(1)+1
      T=TT
      CALL NORM(A,NA,ANAA)
      TMAX= 10.01/ANAA
      K=0
  101 IF (TMAX-T ) 103,104,104
  103 K=K+1
      T=TT/2**K
      IF (K-1000) 101,102,102
  104 SC=T
      CALL SCALE(A,NA,A,NA,T)
      CALL UNITY(B,NB)
      II=2
      N = 35
      CALL ADD(A,NA,B,NB,DUMMY(1),ND)
      CALL EQUATE(A,NA,DUMMY(NDD),ND)
  106 CALL MULT(A,NA,DUMMY(NDD),ND,B,NB)
      S=1.D0/II
      CALL SCALE(B,NB,DUMMY(NDD),ND,S)
      CALL ADD(DUMMY(NDD),ND,DUMMY(1),ND,B,NB)
      CALL EQUATE(B,NB,DUMMY(1),ND)
      N=N-1
      IF (N) 107,107,105
  105 II=II+1
      GO TO 106
  107 IF (K) 109,108,212
  109 CALL LNCNT (1)
      WRITE (6,110)
  110 FORMAT (' ERROR IN ETPHI  K IS NEGATIVE')
  112 IF (K-1) 213,212,212
  213 K=1
  212 DO 111 J=1,K
      T=2*T
      CALL EQUATE(B, NB, DUMMY(1), ND)
      CALL EQUATE(DUMMY(1), ND, DUMMY(NDD), ND)
  111 CALL MULT(DUMMY(NDD),ND,DUMMY(1),ND,B,NB)
      TT=T
  108 CONTINUE
      S=1.D0/SC
      CALL SCALE(A,NA,A,NA,S)
      RETURN
  102 CALL LNCNT (1)
      WRITE (6,119)
  119 FORMAT   (' ERROR IN ETPHI   K=1000')
      CALL ASPERR
      RETURN
  998 CALL LNCNT (1)
      WRITE  (6,50)  NA,KDUM,NDD
   50 FORMAT  (' DIMENSION ERROR IN ETPHI   NA='2I6,   'KDUM='I5,5X,
     1   'KDUM(MIN)='I5)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.AUG$$
FTN AUG$$,N,ISD=N
      SUBROUTINE AUG(F,NF,G,NG,RI,NRI,H, NH,Q,NQ, C,NC,Z,NZ, II)
      DIMENSION F(1),G(1),RI(1),H(1),Q(1),Z(1),C(1)
      DIMENSION NNP1(2),NNP2(2),NNP3(2),NNP4(2),NF(2),NG(2),NRI(2),
     1NH(2),NZ(2),NC(2),NN(2),NQ(2)
      DOUBLE PRECISION F, G, RI,H,Q,C,Z
      IF((NF(1).NE.NF(2)).OR.(NRI(1).NE.NRI(2)).OR.(
     1NQ(1).NE.NQ(2))) GO TO 995
      NNZ=2*NF(1)
      IF( (NG(1).NE.NF(1)).OR.(NG(2).NE.NRI(1)))GO TO 995
      IF(II.EQ.1) GO TO 206
      IF((NH(1).NE.NQ(1)).OR.(NH(2).NE.NF(2))) GO TO 995
  206 TWO = 2
      N = NF(1)
      NSQ = N*N
      NZ(1)=NNZ
      NZ(2)=NNZ
      NP1=1
      NP2 = NP1 + NSQ
      NP3 = NP2+NSQ
      NP4 = NP3 + NSQ
      CALL TRANP(G,NG,Z(NP4),NNP4)
      CALL MULT(RI,NRI,Z(NP4),NNP4,C,NC)
      CALL MULT(G,NG,C,NC,Z(NP4), NNP4)
      IF(II .EQ. 1) GO TO 204
      CALL TRANP(H,NH,Z(NP3), NNP3)
      CALL MULT(Q,NQ,H,NH,Z(1), NNP1)
      CALL MULT(Z(NP3),NNP3,Z(  1),NNP1,Z(NP2),NNP2)
      GO TO 205
  204 CALL EQUATE(Q, NQ, Z(NP2), NQ)
  205 NPAIR=MOD(N,2)
      IF(NPAIR.EQ.0) NPAIR=TWO
      NN(1) = N
      NN(2) = 1
      GO TO (201,202),NPAIR
  201 DO 300 I=1,N,2
      NP2 = N*(N+I-1)+1
      NTH3=TWO*N*(I-1)+N+1
  300 CALL EQUATE(Z(NP2),NN,Z(NTH3),NN)
      DO 302 I=2,N,2
      NP4=N*(3*N+I-1)+1
      NTH2=TWO*N*(N+I-1)+1
  302 CALL EQUATE(Z(NP4),NN,Z(NTH2),NN)
      GO TO (202,203),NPAIR
  202 DO 301 I=2  ,N,2
      NP2 = N*(N+I-1)+1
      NTH3=TWO*N*(I-1)+N+1
  301 CALL EQUATE(Z(NP2),NN,Z(NTH3),NN)
      DO 304 I=1,N,2
      NP4=N*(3*N+I-1)+1
      NTH2=TWO*N*(N+I-1)+1
  304 CALL EQUATE(Z(NP4),NN,Z(NTH2),NN)
      GO TO (203,201),NPAIR
  203 DO 303 I=1,N
      IJ = I+N
      DO 303 J=1,N
      JJ = J+N
      L1=NNZ*(J-1)+I
      L2=NNZ*(IJ-1)+JJ
      L3=N*(J-1)+I
      Z(L1)=-F(L3)
  303 Z(L2)=F(L3)
      GO TO 1000
  995 CALL LNCNT (2)
      WRITE  (6,50)  NF,NG,NRI,NH,NQ
   50 FORMAT  (' DIMENSION ERROR IN AUG',T35,'NF',9X,'NG',9X,'NRI',8X,
     1   'NH',9X,'NQ'/29X,5(3X,2I6))
  999 CALL ASPERR
 1000 RETURN
      END
ERASE SOURCE.RICAT$$
FTN RICAT$$,N,ISD=N
      SUBROUTINE RICAT(PHI,NPHI,C,NC,NCONT,K,NK,PT,NPT,W,KDUM)
      DIMENSION NCONT(3),NPHI(2),NC(2),NK(2),NN(2),NM(2),NPT(2)
      DIMENSION PHI(1),C(1),K(1),PT(1),W(1)
      DOUBLE PRECISION PHI, C, K, PT, SUM, SUMN, AL, W,DET
      TWO=2
      N = NPHI(1)/TWO
      NSQ=N*N
      NSQ4=4*NSQ
      NP1=1
      NP2= NSQ+NP1
      NP3=NSQ+NP2
      NP4= NSQ+NP3
      IF  (KDUM.LT.NSQ4 )  GO TO 600
      IF  (NPHI(2).NE.NPHI(1).OR.NC(2).NE.N.OR.NPT(1).NE.N.OR.NPT(2).
     1  NE.N)   GO TO 600
      NPRINT=NCONT(1)
      NBLOCK=NCONT(2)/NPRINT
      NZ=NCONT(3)
C     REARRANGE PHI MATRIX
      NN(1)=N
      NN(2)=1
      DO 300 I=1,N
      NTH1 =TWO*N*(I-1)+1
      NTH3=NTH1+N
      NW1=N*(I-1)+1
      NW2=NW1+N*N
      CALL EQUATE(PHI(NTH1),NN,W(NW1),NN)
  300 CALL EQUATE(PHI(NTH3),NN,W(NW2),NN)
      NM(1)=TWO*N*N
      NM(2)=1
      CALL EQUATE(W(1),NM,PHI(1),NM)
      DO 301 I=1,N
      NTH2=TWO*N*(N+I-1)+1
      NTH4=NTH2+N
      NW3 = N*(TWO*N+I-1)+1
      NW4= NW3+N*N
      CALL EQUATE(PHI(NTH2),NN,W(NW3),NN)
  301 CALL EQUATE(PHI(NTH4),NN,W(NW4),NN)
      NWW=TWO*N*N+1
      CALL EQUATE(W(NWW),NM,PHI(NWW),NM)
C     MAIN COMPUTATION
C     CALL UNITY(PT,NPT)
      DO 306 I= 1,N
  306 K(I) = 0.
      NTOT=0
      DO 403 I=1,NBLOCK
      DO 104 J=1,NPRINT
      CALL MULT(PHI(NP3), NPT, PT, NPT, W(1), NPT)
      CALL ADD (PHI(1), NPT, W(1), NPT, W(1), NPT)
      CALL INV(W(1), NPT, DET, W(NP2))
      CALL MULT(PHI(NP4), NPT, PT, NPT, W(NP2), NPT)
      CALL ADD(PHI(NP2), NPT, W(NP2), NPT, W(NP2), NPT)
      CALL MULT(W(NP2), NPT, W(1), NPT, PT, NPT)
      SUMN=0.
      SUM=0.
      DO 202 IL=1,N
      ILP=IL+NP3
      NIL=(IL-1)*N+IL
      SUM=SUM+DABS(PT(NIL))
  202 SUMN=SUMN+DABS(PT(NIL)        -W(ILP))
      AL=SUMN/SUM
      DO 201 IL=1,N
      NIL=(IL-1)*N+IL
      ILP=IL+NP3
  201 W(ILP)   =PT(NIL)
  204 DO 104 M=2,N
      N1=M-1
      DO 104 L=1,N1
      L1=N*(L-1)+M
      L2=N*(M-1)+L
      PT(L1)=(PT(L1) + PT(L2))/2.
      PT(L2)=PT(L1)
      IF(AL-.00001) 203,203,104
  104 CONTINUE
      NTOT=I*NPRINT
      GO TO 305
  203 NTOT=NTOT+J
  305 CALL  MULT (C,NC,PT,NPT,K,NK)
  103 GO TO (403,400,401,402), NZ
  400 CALL LNCNT (1)
      WRITE (6, 50) NTOT
   50 FORMAT (10X,I4,14H   ITERATIONS  )
      CALL PRNT  (PT,NPT,'P(T)',1)
      GO TO 403
  401 CALL LNCNT (1)
      WRITE (6, 50) NTOT
      CALL PRNT   (K,NK,'K(T)',1)
      GO TO 403
  402 CALL LNCNT (1)
      WRITE (6, 50) NTOT
      CALL PRNT   (K,NK,'K(T)',1)
      CALL PRNT  (PT,NPT,'P(T)',1)
      IF(AL-.00001) 210,210,403
  403 CONTINUE
C     REARRANGE PHI MATRIX
  210 CALL EQUATE(PHI(1),NM,W(1),NM)
      DO 303 I=1,N
      NTH1 =TWO*N*(I-1)+1
      NTH3=NTH1+N
      NW1=N*(I-1)+1
      NW2=NW1+N*N
      CALL EQUATE(W(NW1),NN,PHI(NTH1),NN)
  303 CALL EQUATE(W(NW2),NN,PHI(NTH3),NN)
      CALL EQUATE(PHI(NWW),NM,W(NWW),NM)
      DO 304 I=1,N
      NTH2=TWO*N*(N+I-1)+1
      NTH4=NTH2+N
      NW3 = N*(TWO*N+I-1)+1
      NW4= NW3+N*N
      CALL EQUATE(W(NW3),NN,PHI(NTH2),NN)
  304 CALL EQUATE(W(NW4),NN,PHI(NTH4),NN)
      RETURN
  600 CALL LNCNT (2)
      WRITE  (6,1600)  NPHI,NC,NPT,KDUM,NSQ4
 1600 FORMAT  (' DIMENSION ERROR IN RICAT',T35,'NPHI',7X,'NC',9X,'NPT'
     1   ,6X,'KDUM',3X,'KDUM(MIN)'/29X,3(3X,2I4),4X,I4,5X,I4)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.SAMPL$$
FTN SAMPL$$,N,ISD=N
      SUBROUTINE SAMPL(PHI,NPHI,H,NH,Q,NQ,R,NR,P,NP,K,NK,NCONT,DUM,KDUM)
      DIMENSION  NPHI(2),NH(2),NQ(2),NR(2),NP(2),NK(2),NCONT(3),NZD(12)
      DOUBLE PRECISION PHI(1),H(1),Q(1),R(1),P(1),K(1),DUM(1)
      DOUBLE PRECISION SUM,SUMN,AL
C     DIMENSION OF DUM MUST BE AT LEAST 6*N*N
C     CHECK FOR CONFORMABLE MATRICES
      N=NPHI(1)
      M=NH(1)
      NK(1)=N
      NK(2)=M
      NSQ=N*N
      ND1=1
      ND2=NSQ+1
      ND6=5*NSQ+1
      ND3=ND2+NSQ
      NSQ6=6*NSQ
      IF (NPHI(2).NE.N.OR.NH(2).NE.N.OR.NQ(1).NE.N.OR.NQ(2).NE.N.OR.NR(1
     1).NE.M.OR.NR(2).NE.M.OR.NP(1).NE.N.OR.NP(2).NE.N.OR.KDUM.LT.NSQ6)
     2GO TO 900
      NFIN=NCONT(2)
      NPRINT=NCONT(1)
      NZ=NCONT(3)
C     START OF MAIN COMPUTATION,  P(0) IS INPUT DATA IN P MATRIX
      KLN=0
      JFN=0
      I=0
      AL=1.
  100 CALL MULT(H,NH,P,NP,DUM,NZD)
C     DUM=H*P
      CALL TRANP (H,NH,DUM(ND2),NZD(3))
C     DUM2=HPRIME
      CALL MULT  (DUM,NZD,DUM(ND2),NZD(3),DUM(ND3),NZD(5))
C     DUM3=H*P*HPRIME
      CALL ADD   (DUM(ND3),NZD(5),R,NR,DUM,NZD)
C     DUM=H*P*HPRIME+R
      CALL PSEUDO (DUM,NZD,DUM(ND2),NZD(3),DUM(ND3),KDUM-3*NSQ)
C     DUM2=INVERSE
      CALL TRANP (H,NH,DUM,NZD)
C     DUM=HPRIME
      CALL MULT  (DUM,NZD,DUM(ND2),NZD(3),DUM(ND3),NZD(5))
      CALL MULT  (P,NP,DUM(ND3),NZD(5),DUM,NZD)
C     DUM=A
  110 CALL MULT  (PHI,NPHI,DUM,NZD,K,NK)
      CALL MULT  (DUM,NZD,H,NH,DUM(ND2),NZD(3))
C     DUM2=A*H
      CALL MULT  (DUM(ND2),NZD(3),P,NP,DUM,NZD)
C     DUM=A*H*P
      CALL SUBT  (P,NP,DUM,NZD,DUM,NZD)
C     DUM= P-A*H*P
      CALL TRANP (PHI,NPHI,DUM(ND2),NZD(3))
      CALL MULT  (DUM,NZD,DUM(ND2),NZD(3),DUM(ND3),NZD(5))
      CALL MULT  (PHI,NPHI,DUM(ND3),NZD(5),DUM,NZD)
C     DUM=PHI(P-A*H*P)PHIPRIME
      CALL ADD   (DUM,NZD,Q,NQ,P,NP)
      DO 120 M=2,N
      N1=M-1
      DO 120 L=1,N1
      L1=N*(L-1)+M
      L2=N*(M-1)+L
      P (L1)=(P (L1)+P (L2))/2.
  120 P (L2)=P (L1)
  130 IF (I.EQ.0)  GO TO 150
      SUM=0.
      SUMN=0.
      DO 140 IL=1,N
      IJ=(IL-1)*N+IL
      SUM=SUM+DABS(P(IJ))
      NDL=ND6+IL
  140 SUMN=SUMN+DABS(P(IJ)-DUM(NDL))
      AL=SUMN/SUM
  150 DO 151 IL=1,N
      IJ=(IL-1)*N+IL
      NDL=ND6+IL
  151 DUM(NDL)   =P(IJ)
      ILST=I
      I=I+1
      IF (AL.LE..00001)  GO TO 300
      IF  (I.GE.NFIN)   GO TO 310
C     INTERMEDIATE PRINT
      IF (I.LT.NPRINT  ) GO TO  100
      NPRINT=NPRINT+NCONT(1)
  152 GO TO  (170,156,155,155),NZ
  155 CALL LNCNT(2)
      WRITE  (6,1152) ILST
 1152 FORMAT  ('0STEP NUMBER=',I4,' IN SAMPL')
      CALL PRNT  (K,NK,4HK(T),1)
      GO TO  (170,170,170,156),NZ
  156 CALL LNCNT(2)
      WRITE  (6,1152) I
  160 CALL PRNT  (P,NP,4HP(T),1)
  170 IF (JFN.EQ.0)  GO TO 100
      RETURN
  300 CALL LNCNT(2)
      WRITE   (6,1300)
 1300 FORMAT   ('0P NO LONGER CHANGING,   EXIT FROM SAMPL')
  310 JFN=1
      GO TO 152
  900 CALL LNCNT (2)
      WRITE (6,1000) NPHI,NH,NQ,NR,NP,KDUM,NSQ6
 1000 FORMAT  ('DIMENSION ERROR IN SAMPL',T35     ,3X,4HNPHI,8X,2HNH,
     19X,2HNQ,9X,2HNR,9X,2HNP,5X,4HKDUM,3X,9HKDUM(MIN)/29X,5(3X,2I4),
     23X,I4,5X,I4)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.TRNSI$$
FTN TRNSI$$, BCD=Y
      SUBROUTINE TRNSI(F,NF,G,NG,J,NJ,R,NR,K,NK,H,NH,X,NX,T,DUMMY,KDUM)
      DOUBLE PRECISION F,G,J,K,H,X,T,Y,R,DUMMY
      DIMENSION F(1),NF(1),G(1),NG(1),J(1),NJ(1),K(1),NK(1),H(1),NH(1),
     1X(1),NX(1),T(1),R(1),NR(1),DUMMY(1)
      DIMENSION NNF(2),NNG(2),NU(2),NNX(2),NY(2)
      DATA STAR/'*'/
      IF(NF(2).NE.NG(1) .OR. NJ(2).NE.NR(1) .OR. NK(2).NE.NX(1) .OR.
     1NJ(1).NE.NK (1) .OR. NR(2).NE.NX(2) .OR.NH(2).NE.NX(1) .OR.
     2NF(2) .NE.NX(1)) GO TO 999
      MAX = NF(1)*(NF(2) +NG(2)+ 1) +NH(1)+NK(1)
      IF (KDUM .LT. MAX) GO TO 910
      I1 = 1
      NSQ =NF(I1) *NF(I1)
      NX4 = NSQ   *4
C
      IF(KDUM .LT. NX4) GO TO 900
C    TRNSI PROGRAM
      N2 = I1 + NSQ
      N3 = N2 + NSQ
      N4 = N3 + NSQ
      L3 = N2 +NF(I1)*NG(2)
      L4 = L3 + NJ(I1)*NR(2)
      L5 = L4 + NH(1)
      L6 = L5 + NJ(I1)*NR(2)
      T1 =T(1)
      N = (T(2) + .5*T1)/T1
      NXR =NX(I1)
      LAST = L6 -I1
      TT = T(4)
      CALL PRNT (F,NF,'  F ',1)
      CALL EAT(F,NF,T1,DUMMY(N3),NNF,DUMMY(N4),NNF,DUMMY(I1),KDUM)
  100 FORMAT(1H0 1P8E16.7)
      CALL  PRNT(DUMMY(N3), NF, 'EAT ', 1)
      CALL EQUATE(DUMMY(N3),NF,DUMMY(I1),NNF)
      CALL  PRNT(DUMMY(N4), NF, 'INT ', 1)
      CALL MULT(DUMMY(N4), NF,G,NG,DUMMY(N2),NNG)
      CALL MULT(J,NJ,R,NR,DUMMY(L3),NU)
      CALL LNCNT (100)
      CALL LNCNT (3)
      WRITE(6, 50)
   50 FORMAT(1H0 'TRANSIENT RESPONSE, * INDICATES CONTROL CHANGES')
      WRITE(6, 51) NXR,    NH(I1),NK(I1)
   51 FORMAT(1H0 4X,'TIME    FIRST',I3,'  ELEMENTS CONTAIN X, NEXT',I3,'
     1  ELEMENTS CONTAIN Y = HX, LAST',I3,'  ELEMENTS CONTAIN U =JR -KX'
     2)
  201 NQ=1
      IU = I1
      CALL MULT(K,NK, X,NX, DUMMY (L5), NU)
      CALL SUBT(DUMMY(L3), NU,DUMMY(L5), NU, DUMMY(L5),NU)
  203 CALL MULT(H,NH,X,NX,DUMMY(L4 ),NY)
      IF (IU .NE. I1) GO TO 205
      WRITE(6, 101) STAR,TT,(X(IP), IP=1,NXR),(DUMMY(IP),IP=L4,LAST)
  101 FORMAT(1H0 A1,F8.2,1P7D15.7/(10X,1P7D15.7))
      GO TO 206
  205 WRITE(6, 102)      TT,(X(IP), IP=1,NXR),(DUMMY(IP),IP=L4,LAST)
  102 FORMAT(1H0 1X,F8.2,1P7D15.7/(10X,1P7D15.7))
  206 IU = I1  +IU
      IF (TT .GE. DABS( T(3)))  RETURN
      TT = TT + T1
  202 CALL MULT(DUMMY(I1), NNF, X, NX, DUMMY(L6 ),NNX)
      CALL MULT(DUMMY(N2),NG, DUMMY(L5),  NU, X, NNX)
      CALL ADD(X,NX,DUMMY(L6 ),NNX, X,NNX)
      IF(NQ .GE. N) GO TO 201
      NQ = NQ + 1
      GO TO 203
C
C    DIMENSION ERROR DIAGNOSTIC
  999 WRITE(6, 990)
  990 FORMAT(1H0 'DIMENSION ERROR IN TRNSI'/25X,'COL SIZE OF 1ST MATRIX
     1     ROW SIZE OF 2ND MATRIX')
      WRITE(6,991) NF(2), NG(1)
      WRITE(6,992) NJ(2),NR(1)
      WRITE(6,993) NK(2),NX(1)
      WRITE(6,995) NH(2),NX(1)
      WRITE(6,996) NF(2), NX(1)
      WRITE(6,994) NJ(1),NK(1),NR(2),NX(2)
  991 FORMAT(1H0 'INTEGRAL (EXP(F.T)) G ' I13,20X, I8)
  992 FORMAT(1H0 'J R' 17X,I15,20X,I8)
  993 FORMAT(1H0 'K X' 17X,I15,20X,I8)
  994 FORMAT(1H0 'JR -KX     ROW SIZE OF J IS'I5,3X,'OF K IS',I5,3X,'COL
     1 SIZE OF R IS' I5,3X, 'OF X IS' I5)
  995 FORMAT(1H0 'H X' 17X,I15,20X,I8)
  996 FORMAT(1H0 'EXP(F.T) X' 10X,I15,20X,I8)
      GO TO 1000
  900 WRITE(6, 52) NX4,KDUM
   52 FORMAT(1H0 'DUMMY MUST BE DIMENSIONED AT LEAST' I4, ' X 1' 'BUT IS
     1 DIMENSIONED ONLY' I4,' X 1')
      GO TO 1000
  910 WRITE(6, 52) MAX,KDUM
 1000  CALL  ASPERR
      RETURN
      END
ERASE SOURCE.PSEUDO$$
FTN PSEUDO$$,N,ISD=N
      SUBROUTINE PSEUDO(A,NA,B,NB,DUM,KDUM)
      DIMENSION A(1),B(1),DUM(1),NA(2),NB(2),IP(4),DEP(3)
      DOUBLE PRECISION A,B,DUM,DEP
      NNW=3*NA(1)*NA(2)
      DO 10 I=1,2
      DEP(I)=0.0
   10 IP(I)=0
   20 IP(3)=NA(1)
      IP(4)=NA(2)
      NB(1)=NA(2)
      NB(2)=NA(1)
      IF  (KDUM.LT.NNW)  GO TO 999
      NEE=NA(1)*NA(2) + 1
      ND=2*NEE - 1
C  IF A IS (1,1) MATRIX ROUTINE INVERTS A(1) AND PUTS IT IN B(1)
      IF  (NA(1).EQ.1.AND.NA(2).EQ.1)   GO TO 600
      CALL PSEU(A,B,DUM,DUM(NEE),DEP,IP,DUM(ND))
      GO TO 1000
  600 B(1)=1./A(1)
      GO TO 1000
  999 WRITE(6,500) KDUM,NNW
  500 FORMAT  (' DIMENSION ERROR IN PSEUDO   KDUM='I5,3X,'KDUM(MIN)='I5)
 1000 RETURN
      END
ERASE SOURCE.DECGEN$$
FTN DECGEN$$,N,ISD=N
C  SUBROUTINE DECGEN     DECOMPOSE A REAL    GENERAL MATRIX
      SUBROUTINE DECGEN (R,NR,G,NG,H,NH,DUM,KDUM)
      DIMENSION  NR(2),NH(2),NG(2),ND(2),NA(2),KP(4)
      COMMON  /MAX/MAXRC
      DOUBLE PRECISION  R(1),G(1),H(1),DUM(1),DCM(3)
      NJ=NR(1)*NR(2)+1
      CALL TRANP (R,NR,DUM(NJ),ND)
      IF(NR(1).GT.NR(2))GO TO 10
      CALL MULT (R,NR,DUM(NJ),ND,DUM,NA)
      ICS=1
      GO TO 30
   10 CALL MULT (DUM(NJ),ND,R,NR,DUM,NA)
      ICS=2
      GO TO 30
C  ENTRY DECSYM         DECOMPOSE A REAL SYMETRIC MATRIX
      ENTRY DECSYM (R,NR,G,NG,H,NH,DUM,KDUM)
      ICS=0
      IF  (NR(1).NE.NR(2))  GO TO 601
      CALL EQUATE (R,NR,DUM,NA)
C       DUMMY STORAGE  MUST BE AT LEAST 7*NA(1)**2
   30 M=NA(1)
      MM7=7*M*M
      IF  (KDUM.LT.MM7  )  GO TO 601
      KP(1)=0
      KP(2)=0
      KP(3)=M
      KP(4)=M
      DCM(1)=0.
      DCM(2)=0.
      MS=M*M
      N2=MS+1
      N3=N2+MS
      N4=N3+MS
      N5=N4+MS
      N6=N5+MS
      N7=N6+MS
      CALL DECOM (DUM,DUM(N7),DUM(N4),DUM(N5),DUM(N6),DCM,KP,DUM(N2))
      CALL TRANP  (DUM(N3),NA,DUM(N6),ND)
      IF (ICS.EQ.2)  GO TO  100
      CALL MULT (DUM(N2),NA,DUM(N6),ND,H,NH)
      NH(2)=KP(2)
      IF (ICS.EQ.1)  GO TO 60
      CALL TRANP (H,NH,G,NG)
      GO TO 150
   60 CALL MULT  (DUM(N5),NA,H,NH,G,NG)
      CALL TRANP (G,NG,DUM(N6),ND)
      CALL MULT  (DUM(N6),ND,R,NR,G,NG)
      GO TO 150
  100 CALL MULT (DUM(N2),NA,DUM(N6),ND,H,NH)
      NG(2)=KP(2)
      CALL TRANP  (H,NH,G,NG)
      CALL MULT  (G,NG,DUM(N5),NA,H,NH)
      CALL TRANP  (H,NH,DUM(N6),ND)
      CALL MULT  (R,NR,DUM(N6),ND,H,NH)
  150 DUM(N6)=KP(2)
      RETURN
  601 CALL LNCNT (1)
      WRITE  (6,1601) NR,KDUM,MM7
 1601 FORMAT(' DIMENSION ERROR IN DECSYM (DECGEN)   NR='2I6,3X,
     1    'KDUM='I4,3X,'KDUM(MIN)='I4)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.READ1$$
FTN READ1$$,N,ISD=N
      SUBROUTINE READ1 (A,NA,NZ,NAM)
      COMMON  /MAX/MAXRC
      DIMENSION  A(1  ),NA(2),NZ(2)
      DOUBLE PRECISION A
      IF  (NZ(1).EQ.0)  GO TO 410
      NR=NZ(1)
      NC=NZ(2)
      NLST=NR*NC
      IF(NLST .GT. MAXRC .OR. NLST .LT. 1.OR.NR.LT.1)  GO TO 16
      DO 400 I = 1, NR
  400 READ (5,101) (A(  J), J = I,NLST,NR)
      NA(1)=NR
      NA(2)=NC
  410 CALL  PRNT (A,NA,NAM,1)
  101 FORMAT (8F10.2)
      RETURN
   16 CALL LNCNT(1)
      WRITE  (6,916)  NAM,NR,NC
  916 FORMAT  (' ERROR IN READ1   MATRIX 'A4,' HAS NA=',2I6)
      CALL ASPERR
      RETURN
      END
ERASE SOURCE.ASPERR$$
FTN ASPERR$$,N,ISD=N
      SUBROUTINE ASPERR
      DATA I /10/
C    ERRTRA IS THE 360/67 TRACE ROUTINE     TRACE IS FOR TSS
C     CALL TRACE
      CALL ERRTRA
C        ERRTRA IS PART OF THE IBM EXTENDED ERROR HANDLING FACILITY
C        TO PROVIDE AN ERROR WALKBACK
C        IF THE PROGRAM IS USED UNDER TSS/360 EITHER OMIT THE
C        CALL ERRTRA COMPLETELY OR MAKE IT A COMMENT  - AND USE
C        CALL TRACE
      I=I-1
      IF (I.GT.0)  RETURN
      I=10
      WRITE (6,100)
  100 FORMAT (' TOO MANY ERRORS.   EXIT CALLED')
      CALL EXIT
      RETURN
      END
ERASE SOURCE.BLKDTA$$
FTN BLKDTA$$,N,ISD=N
      BLOCK DATA
      COMMON  /FORM/NEPR,FMT1(6),FMT2(6)
      COMMON/LINES/NLP,LIN,TITLE(23)
      COMMON  /MAX/MAXRC
      DATA  MAXRC/6400/
C- NOTE NLP NO. LINES/PAGE VARIES WITH THE INSTALLATION.
      DATA LIN,NLP/1,45/
      DATA  NEPR,FMT1     /7,'(1P7D16.7)'/
      DATA  FMT2/'(3X,1P7D16.7)'/
      DATA  TITLE  /19*'    ','VASP PROGRAM    '/
      END
ERASE SOURCE.PSEU$$
FTN PSEU$$,N,ISD=N
      SUBROUTINE  PSEU(A,B,C,EE,DEP,IP,  D)
C SUBR PSEU GENERALIZED INV BY ANDR.  ALG. J ANDREWS, INF SYS CO. APR 69
C- SUBR TO COMPUTE PSEUDO-INVERSE OF GENERAL MATRIX, RETURN FINAL PIVOT
C... NOTE IMPLIT STATEMENTS MUST BE -FIRST- CAN BE REPLACED BY TYPE
      IMPLICIT REAL*8 (D), INTEGER*2 (Q)
      COMMON  /MAX/MAXRC
C DOUBLE PRECISION IS THE ONLY THING ESSENTIAL.
      INTEGER*2 M
      DOUBLE PRECISION A,B,C,EE, D
      DIMENSION A(400),B(400),C(400),EE(400),  D(2000),
     1 KRV(4),
     2 DEP(3), DPR(2), IP(4), JP(5)
      DATA ICC, DFZO /Z40000000, 0.D0 /
      EQUIVALENCE(DDI,FDI,IDD),(DMX,FMX)
      EQUIVALENCE(DDI,DSUM),(DFZO,FZRO,IZ,QZ),(QLL,QR1), (KRV(1),KRC),
     1 (KRV(2),KRC2),(KRV(3),KRC3),(KRV(4),KRC4)
      QPS = 1
      GO TO 1000
       ENTRY    PSEU  P(A,B,C,EE,DEP,IP,  D)
      QPS = IZ
      IP(4) = IP(3)
1000  CONTINUE
      DP1 = DEP(1)
      EF2 = SNGL(DEP(2))
C- SET DEFAULT VALUES OF TOLERANCES
      IF(DEP(1) .EQ. DFZO) DP1 = 2.D-6
      IF(EF2 .EQ. FZRO) EF2 = 1.0
      NCA = IP(4)
C NUMBER OF ROWS OF ORIGNAL INPUT MATRIX
      QR = IP(3)
C- SET SW FOR =0, DO ALL STEPS, NOT=0, THEN WANT RANK ONLY.
      QNT =  QR*NCA
C- TEST DIMENSIONS INPUT FOR REASONABLENESS.
      IF(QNT .LT. 2 .OR. QNT .GT. MAXRC.OR.QR.LT.1)  GO TO 691
C- IF DIMENSIONS ABSURD, PSEU ERR EXIT 1.
      QDCM = IP(1)
      QITR = QDCM
      IF(QDCM .LT. QZ) QITR = QDCM +1
      NR = QR
      QIT = IZ
C- TEST TO SEE IF SYMMETRIZATION IS NEEDED.
      IF(QPS ) 16, 150, 16
C- TEST TO FIND SMALLER DIMENSION OF MATRIX.
16    IF(QR - NCA) 18,18,19
19    NR = NCA
      QX = 1
      QR2 = QR
      QLL = QR
      QTP = IZ
      GO TO 170
18    CONTINUE
      QX = NR
      QR2 = 1
      QLL = NCA
      QTP  = 1
170   CONTINUE
C- SET ROW-COLUMN LIMIT TO APPROPRIATE CASE, EITHER ROW OR COLM DIMENS.
      DO 181 I =1, NR
      DO 181 K =1, NR
      DSUM = DFZO
      DO 183 J = 1, QLL
      QN = QX*J - QR
      L = QN + QR2*I
      M = QN + QR2*K
183   DSUM = DSUM + A(L)*A(M)
C- SUM OF A(I,LL) * A(K,LL),, LL RUNS 1 TO ROW OR TO COLM LIMIT
C. B = A*A TRANS HAS COLM LIMIT, B = A TRANS * A HAS ROW LIMIT.
      L = (K-1)*NR + I
181   B(L) = DSUM
      GO TO 188
C- HERE MOVE A TO B. A IS ALREADY POSITIVE DEFINITE.
150   DO 151 L = 1, QNT
151   B(L) = A(L)
C-. FORCE SYMMETRIZATION OF B, TO COMPENSATE FOR ROUND-OFF,MULTIPLIC.
188   DO 189 I =1, NR
      DO 189 K =1, NR
C,  B(I,K) = B(K,I) = 1/2 (B(I,K) + B(K,I) )
      L = I + (K-1)*NR
      M = K + (I-1)*NR
      DSUM = (B(L) + B(M) ) * 0.5D0
      B(L) = DSUM
189   B(M) = DSUM
C HERE SET UP CALL -INITIAL- OF ANDRA. ONLY COMES HERE ONCE PER MATRIX.
      QNT = NR*NR
      KRC = QNT
      KRC2 = QNT + KRC
      KRC3 = QNT + KRC2
      KRC4 = QNT + KRC3
C-* OMIT SAVING  OF B, IF RANK ONLY AND NO ITERATION
      IF(IP(2) .NE. IZ .AND. QITR .EQ. IZ) GO TO 200
      DO 1891 I =1, QNT
1891  D(I) = B(I)
200   CONTINUE
C- SEARCH DIAGONAL OF INPUT FOR LARGEST ELEMENT. USE TO DEFINE FL. PT.
      QR1 = NR + 1
      L = 1
      DMX = DFZO
      M = IZ
      DO 23 I = 1, NR
      DDI = DABS( B(L) )
      IF(FMX .GE. FDI) GO TO 23
      M = L
      FMX = FDI
23    L = QR1 + L
      IF(M .EQ. IZ) GO TO 692
C- SET TOLERANCE FOR ANDRA  LIMIT OF SIZE OF DIAGONAL.
C TOLERANCE OF ZERO IN ANDRA CALL.
      DPR(1) = DABS(DP1* B(M))
C - ASK FOR ALL ROWS, DONE IN 1 CALL
      JP(1) = IZ
C- JP2 FIRST TIME INITIALIZATION FOR ANDRA
      JP(2) = IZ
      IF(QIT .NE. QZ) GO TO 561
      JP(4) = NR
      M = IZ
      SOLD = - EF2
C --  HAVE FINISHED PRELIM. PART
C  INITIALIZATION FOR ANDRA (DIAGONALIZATION) NOW COMPLETED.
CALL ANDRA  -TO DIAGONALIZE SYMMETRIC MATRIX.
C CALL ANDRA REDUCES ROWS BY MODIFIED GAUSS METHOD, USING SQRT(PIVOT).
30    CONTINUE
      IF(QITR .EQ. QZ) GO TO 32
C- SAVE OLD VALUES  IN CASE  PIVOT IS REJECTED, UNDER ITERATION OPT.
      DO 31 L =1, QNT
      J = KRC + L
      K = KRC2 + L
      D(J) = B(L)
31    D(K) = C(L)
 32    CALL ANDRA(B,C,DPR,JP)
      JP(1) = QITR
      IR = JP(3)
C- CHECK COMPLETION- IS MATRIX ALL DONE IS MATRIX INVERTIBLE..
      IF(QITR .EQ. IZ .OR. IR .EQ. NR .AND. QIT .EQ. IZ) GO TO 700
CHECK IF ITERATING WITH RHO TEST OR NOT
C* Q1IT IF NO ITERATION  OR NO NEW PIVOT FOUND
C- OMIT ITERATION CALCS. IF NO NEW PIVOT. DECREASE TOLERANCE
      IF(JP(5) .EQ. M) GO TO 220
C COMPUTE RHO FOR ESTIMATING THRESHHOLD TO STOP  SS IS RHO
      SS = (BDNRM(NR,C,EE,D,KRV) + BDNRM(-NR,C,EE,D,KRV) ) *EF2  ** IR
C WHY ONLY SNGLE PREC./THIS IS ONLY A ROUGH TEST TO STOP ITERATION.
C THAT-S WHY. SIMILARLY, OTHER USES OF SINGLE PREC.
      IF(SOLD .LT. SS .AND. SOLD .GT. FZRO) GO TO 650
C- IF SUBSTANTIAL IMPROVEMENT TRY AGAIN,
C, OTHERWISE QUIT, RETURN THE A PSEUDO INVERSE, EVEN IF OFF.
220   CONTINUE
      QIT = QIT +1
      SOLD = SS
C/ SAVE PREVIOUS ROW  IN WHICH A PIVOT WAS FOUND
      M = JP(5)
      IF(QIT .EQ. NR)  GO TO 700
C- PUT IN SMALLER TOLERANCE IN CASE DIAGONAL TOO SMALL OTHERWISE.
      DPR(1) = DPR(2) * 2.D-5
C- TRY TO REDUCE 1 MORE ROW.
      IF(IR - NR) 30, 700,  606
650   CONTINUE
C* RESTORE B AND C TO THEIR PREVIOUS VALUES. THE LAST PIVOT HAS BEEN
CREJECTED (BACK-TRACK), WHILE ITERATING.
      JP(3) = JP(3)  -1
      DO 653 I =1, QNT
      J = KRC + I
      K = KRC2 + I
      B(I) = D(J)
653   C(I) = D(K)
700   CONTINUE
      IR = JP(3)
      M = IZ
C- HERE WISH TO REPLACE MARKERS IN DIAGONAL WITH LEGITIMATE 1.D0
      L = 1
      DO 704 I = 1, NR
      DDI = B(L)
      IF(IDD) 701, 702, 701
701   IF(IDD .NE. ICC) GO TO 7101
      B(L) = 1.D0
      GO TO 704
C AT 7101 FORCE SMALL TRASH TO ZERO.
7101  B(L) = DFZO
702   M = 1
704   L = QR1 + L
C - IF ALREADY TRIED ANOTHER REDUCTION, TO GET MATRIX IN -UPPER- DIAG.
CQR  OMIT PART OF CALCULATIONS IF ONLY RANK IS DESIRED.
      IF(IP(2) .NE. IZ) GO TO 877
C9 QDCM SUPPRESSES LAST PHASE IF DECOM WAS CALLER..
      IF(M .LT. 1 .OR. QDCM .LT. QZ) GO TO 80
C BELOW HAVE SING. MATRIX THAT NEEDS FURTHER WORK.
C- HAVE MATRIX DIAGONALIZED WITH 1S, 0S INTERSPERSED  (A IS SINGULAR)
C- RE-DO TO GET PS-INV THAT MOVES ALL 1S OF DIAGONAL TO UPPER LEFT DIAG.
C.TO COMPUTE U MATRX AS IN ASP, FOR TRANSFORMING ORIG B IN SINGULAR CASE
      L =1
      DO 527 I =1,NR
      DO 525 J =1, NR
      K = (J-1)*NR + I
      IF(B(L) ) 521,522,521
522   C(K) = - C(K)
      C(L) = DFZO
      GO TO 525
521   C(K) = DFZO
      C(L) = 1.D0
525   CONTINUE
527   L = QR1 + L
C-SAVE RANK SO FAR, SHOULD BE SAME SIZE AFTER RE-INVERSION
      QR2 = IR
      DO 54 I = 1, NR
      DO 54 K =1, NR
      DSUM = DFZO
      QN = (K-1)*NR
      DO 53 J =1, NR
      M = (I-1)*NR + J
      L = QN + J
53    DSUM = C(M)*D(L) + DSUM
      L = QN + I
      EE(L) = DSUM
54    CONTINUE
      DO 56 I =1, NR
      DO 56 K =1, NR
      DSUM = DFZO
      QN = (K-1)*NR
      DO 55 J =1, NR
      L = (J-1)*NR + I
      M = QN + J
55    DSUM = EE(L)*C(M) + DSUM
      L = QN + I
      B(L) = DSUM
56    CONTINUE
C, SET UP FOR SECONDARY ANDRA CALL  NO ITERATION   JP4 = NR
      QIT =1
C GO FIND LARGEST DIAG. ELEMENT AGAIN
      GO TO 200
561   JP(3) = IZ
       CALL ANDRA(B,EE,DPR,JP)
      IR = JP(3)
C- TEST FOR A CHANGE IN RANK ... ERROR
      IF(QR2 - IR) 693, 568, 694
568    CALL  T TRM(NR,EE,D)
C- TRANSFORM C SHARP IN D.. BS = ((U)* D *(U TRP) )
      DO 58 I =1, NR
      DO 58 K =1, NR
      DSUM = DFZO
      QN = (K-1)*NR
      DO 57 J =1, NR
      M = (J-1)*NR + I
      L = QN + J
57    DSUM = C(M)*D(L) + DSUM
      L = QN + I
      B(L) = DSUM
58    CONTINUE
      DO 60 I =1,NR
      DO 60 K = 1, NR
      DSUM = DFZO
      DO 59 J =1, NR
      QN = (J-1)*NR
      M = QN + K
      L = QN + I
59    DSUM = B(L)*C(M) + DSUM
      L = (K -1)*NR + I
      EE(L) = DSUM
60    CONTINUE
C- NOW RE-ENTER MAIN SEQUENCE WITH PS-INV. IN EE.
      GO TO 808
C GO FIX UP B PSUEDO-INVERSE. PRESUMABLY HAVE DIAGONALIZED
C HAVE DIAGONALIZED WITH ALL 1S IN UPPER LEFT
C- HERE WE HAVE FINISHED DIAGONALIZ. WANT TO GET PSUEDO INV. IN B.
870   IF(QDCM .LT. QZ) GO TO 877
C NEED TO SAVE DIAGONALIZED B FOR USE BY DECOM CALL (QDCM NEG. FLAG)
      DO 871 I = 1,QNT
C- A WAS SYMMETRIC.  JUST MOVE EE TO B  RETURN FROM PSEUP ENTRY
871   B(I) = EE(I)
      GO TO 877
80    CONTINUE
C NOW FORM (T TRP) * T = APPROX B SHRP PSUEDINV IN MATRX EE
      CALL T TR M(NR,C,EE)
808   IF(QPS .EQ. QZ) GO TO 870
      IF(QTP) 819, 819,818
C HERE B = (A TRANS)*E = (A TRP)*(A*A TRP)-SHRP  NRA .LE. NCA
818   DO 8181 I = 1,NCA
      DO 8181 J = 1,NR
      DSUM = DFZO
      QN = (J -1)*NR
      DO 8182 K = 1,NR
      L = (I -1)*QR + K
      M = QN + K
8182  DSUM = DSUM + A(L)*EE(M)
      L = (J-1)*NCA + I
8181  B(L) = DSUM
      GO TO 877
819   DO 8191 I = 1,NR
      DO 8191 J =1, QR
      DSUM = DFZO
      DO 8192 K = 1,NR
      L = (K -1)*QR + J
      M = (K -1)*NCA + I
8192  DSUM = DSUM + A(L)*EE(M)
C- NOTE NCA IS USED, BECAUSE A-SHARP IS TRANSPOSED IN DIMENSIONS
      L = (J -1)*NCA + I
C HERE B = EE (A TRANS) = (A TRP*A)-SHRP * (A TRANS)  NRA .GT. NCA
8191  B(L) = DSUM
C- HERE GET READY TO RETURN
877   CONTINUE
C- MOVE RANK TO RETURN PARAMETER
      IP(2) = IR
      DEP(3) = DPR(2)
C. ABOVE RETURN FINAL PIVOT FROM ANDRA ALG. DIAGONALIZATION
      RETURN
  691 CALL LNCNT(1)
      WRITE (6,1691) QR,NCA
 1691 FORMAT (' DIMENSION ERROR IN PSEU NA='2I6)
      GO TO 1700
  692 CALL LNCNT(1)
      WRITE (6,1692)
 1692 FORMAT (' ERROR IN PSEU - DIAGONAL ELEMENTS OF MATRIX=0')
      GO TO 1700
  693 CALL LNCNT(1)
      WRITE (6,1693)
 1693 FORMAT (' ERROR IN PSEU   RANK HAS DECREASED   COMPUTATION ENDED')
      GO TO 1700
  694 CALL LNCNT(1)
      WRITE (6,1694)
 1694 FORMAT (' ERROR IN PSEU-RANK HAS INCREASED-COMPUTATION CONTINUES')
      CALL ASPERR
      GO TO 568
  606 CALL LNCNT(1)
      WRITE  (6,1606)
 1606 FORMAT  (' ERROR IN PSEU  RANK IS GREATER THAN MATRIX SIZE')
 1700 CALL ASPERR
      RETURN
      END
ERASE SOURCE.BDNRM$$
FTN BDNRM$$,N,ISD=N
      FUNCTION BDNRM(NR,CT,EE,D,KRV)
      INTEGER*2 QF
      DOUBLE PRECISION CT,EE,D, AN,BR, DFZO,DSUM
      DIMENSION CT(400),EE(400), D(2000), NV(2), KRV(4)
C- D HOLDS 5 MATRICES. THE FIRST AND THE LAST 2 ARE USED HERE
      DIMENSION PPP(2)
      EQUIVALENCE(AN,FN), (BR,FR)
      DATA DFZO /0.D0/
C, EQUIVALENCES BELOW JUST TO SAVE STORAGE
      EQUIVALENCE(DFZO,IZ),(AN,DSUM),(BR,PPP(1),I),(PPP(2),K),
     1 (NV(1),L),(NV(2),M),(IR,NL)
C TEST,, IF NR NEG.,THEN TRANSPOSE ROLES OF D      AND (CT TRANS  *CT)
      QF = NR
      KD3 = KRV(3)
      KD4 = KRV(4)
      IF(NR) 10,10, 20
      ENTRY TTRM(NR,CT,EE)
C TO DO T TR * T ONLY ENTRY TTRM
      QF = IZ
      GO TO 20
10    NR = -NR
20    IR = NR
      DO 30 I = 1, IR
      LL = (I-1)*IR
      DO 30 K = 1, IR
      DSUM = DFZO
      KK = (K -1)*IR
      DO 29 J = 1, IR
      L = J + LL
      M = J + KK
29    DSUM = DSUM + CT(L)*CT(M)
C ABOVE FORMAING T TRANSPOSE TIMES T. WHICH  IS APPROX . OF B SHARP
      L = I + KK
      IF(QF) 31, 39, 32
31    KK = KD3 + L
      D(KK) = D(L)
      EE(L) = DSUM
      GO TO 30
C-39 COMPUTE T TRANSPOSE * T ONLY.. PROVIDES INVERSE B SHARP
39    EE(L) = DSUM
      GO TO 30
32    EE(L) = D(L)
      KK = KD3 + L
      D(KK) = DSUM
30    CONTINUE
      IF(QF) 41,66,41
41    NV(1) = IR
      NV(2) = IR
C- LET P = 1ST MATRIX = EE,, Q = 2D = D(KD3+1)
CCC ROLES OF P AND Q ARE GIVEN TO 2 MATRICES AT 31, SWITCHED AT 32.
COMPUTE D(K4) = Q*P, EE=D(K4)*Q, EE = EE -D(K3) = Q*P*Q - Q
C-FUNCTION IS NRM(Q*P*Q -Q)/NRM(Q)  RESULT = A SCALAR
      CALL MULT(D(KD3+1),NV,EE,NV,D(KD4+1),NV)
      NL = IR*IR
      CALL MULT(D(KD4 +1),NV,D(KD3+1), NV,EE,NV)
      DO 8 I = 1, NL
      KK = KD3 + I
8     EE(I) = EE(I) - D(KK)
      CALL NORM(EE,NV,AN)
      CALL NORM(D(KD3+1),NV,BR)
CQUOTIENT NEARS 0.0 AS BSHRP APPROACHES THAT FITTING 2 MOOR-PENRSE AXIOM
      BDNRM = FN / FR
9     RETURN
   66 BDNRM=FN
C66 IS A DUMMY REALLY WANT  MATRX MULT. ONLY.
C             OS/360 USERS REPLACE BDNRM=FN WITH TTRM=FN
      GO TO 9
C  SIDE COMPUTATIONS   J W ANDREWS  INF. SYSTEMS CO. MAY 1969
      END
ERASE SOURCE.ANDRA$$
FTN ANDRA$$,N,ISD=N
      SUBROUTINE ANDRA(B,T,DPR,JP)
C- SUBROUTINE ANDRA DIAGONALIZES POS.DEF.SYMM. J ANDREWS  I. S. CO.
C - SUBR ANDRA CALLED BY PSEU J W ANDREWS, INF. SYSTEMS CO. APRIL 1969
      IMPLICIT REAL*8 (D), INTEGER*2 (Q)
      DOUBLE PRECISION B, T
      DIMENSION B(400), T(400), DPR(2), JP(5)
      EQUIVALENCE(DDI,FDI,IDD),(DCC,ICC),(DMX,FMX),(DRS,IIS)
      EQUIVALENCE (DFZO,FZRO,IZ)
      DATA ICC, DFZO /Z40000000, 0.D0/
C-DPR1 IS MAGNITUDE THAT IS CONSIDERED ZRO PIVOT MUST BE NO SMALER.
C- DPR(2) IS TO RETURN FINAL PIVOT, SO THAT USER MAY TEST SMALLNESS.
CC- ANDRA CAN BE USED ALL BY ITSELF TO GET INV., RANK OF POS SYMM.
C NOTE THAT DSQRT HAS TO BE TAKEN OF PIVOTS ALONG THE DIAGONAL.
C- NOTE I AM DELIBERATELY ALLOWING SOME PARAMETERS TO CHANGE ON SUBSE-
C-QUENT CALL  DPR(1) CHANGES PIVOT SIZE  A ROUGH TOLERANCE FOR ZRO.
      EF = SNGL(DPR(1) )
C- TEST- IS THIS AN INTIALIZATION CALL/
      IF(JP(2)) 2, 1, 2
C INTIALIZE- FORM IDENTITY MATRIX
1     QS = JP(4)
      QNT = QS*QS
      IF(QS .LT. 1 .OR. QNT .GT. 6400) GO TO 691
      DO 18 I = 1, QNT
18    T(I) = DFZO
      L = 1
      QR1 = QS +1
      DO 1810 I = 1, QS
      T(L) = 1.D0
1810  L = QR1 + L
      DPR(2) = DFZO
C- SET RANK TO ZRO. TRIAL PIVOT VALUE TO ZERO.
      QKR = IZ
C  SET PIVOT CHOICE ITERATION AT 0  ALLOWANCE OF NO. ROWS+1  ITER.
      QITR = IZ
2     CONTINUE
200    CONTINUE
C- ZERO OUT MAX DIAG. AND CT DIAG. TEMPORARY VARIABLES
      FMX = FZRO
      I = IZ
      M = IZ
C- BELOW SEE IF ALL DIAG ELEMENTS TESTED YET
      L = 1 - QR1
30    IF(I .EQ. QS) GO TO 40
      I = I +1
      L = QR1 + L
      DDI = B(L)
C- GET CURRENT DIAG. ELEMENT FOR INTEGER, SINGLE PREC. TEST
C- UPDATE L TO GET -NEXT- DIAG. ELEMENT
C-BELOW TEST FOR DIG. ELEMENT @ALREADY@ REDUCED TO 1.(CODEMARKED),ICC
      IF(IDD .EQ. ICC) GO TO 30
      IF(FDI - FMX) 30,30, 32
C- TEST FOR NEGLIGIBLE FL. PT. QTY.-TREAT THESE, AND NEG., AS ZEROS.
32    IF(FDI .LT. EF) GO TO 30
C- SET NEW MAX, 2BLE PREC., SAVE BEST ROW FOR PIVOT  @QMR@
      QMR = I
      DMX = DDI
      M = L
      GO TO 30
40    CONTINUE
400   IF(M .EQ. IZ) GO TO 505
      DRS = 1. D0 /DSQRT(DMX)
C SET INDEX OF FIRST ROW, QMR (PIVOT) COLUMN
      K = (QMR -1)*QS  +1
      L = QMR
      DO 41 I = 1, QS
      DDM = B(L)*DRS
      T(L) = T(L)*DRS
      B(L) = DDM
C--SYMMETRICALLY, FORCE COLUMN TO SAME VALUE IN B ONLY
      B(K) = DDM
      K = K +1
41    L = QS + L
C FORCE PIVOT ELEMENT TO EXACT VALUE OF UNITY
      B(M) = 1.D0
C- NOW REDUCE ALL OTHER ROWS OF B, T, ELIMINATING COLUMN OF PIVOT VARIAB
      DO 460 I = 1, QS
C- TEST FOR PIVOTAL ROW. OTHER ROWS
      IF(I .EQ. QMR) GO TO 460
COEFF, TO BE ZEROED  CAN NOT BE PREVIOUS PIVOT.
      J = I - QS
      K = QS*I + J
      DRS = B(K)
C  BELOW TEST FOR A ROW ALREADY REDUCED, TO SKIP
      IF(IIS .EQ. ICC) GO TO 460
C- GET COEFF IN PIVOT COLUMN TO BE ELIMINATED
      K = QMR*QS + J
      DMM = - B(K)
      L = QMR
      K = I
      DO 47 J =1, QS
C- L IS ROW USED TO REDUCE, WITH PIVOT.
C  K IS CURRENT ROW THAT PIVOT GETS ELIMINATED FROM.
      B(K) = B(K)  + B(L)*DMM
      T(K) = T(K) + T(L)*DMM
      L = QS + L
47    K = QS + K
460   CONTINUE
      L = QMR
      DO 461 I =1, QS
C  FORCE MOST OF PIVOT ROW TO ZERO. COMPLETES REDUCTION WITH 1 PIVOT/
      B(L) = DFZO
461   L = QS + L
C FORCE PIVOT TO @CODE@ FOR ONE  ..
      B(M) = DCC
C- SIGNAL NO LONGER FIRST TIME CALLED.
      JP(2) = 1
C--   UPDATE EFFECTIVE RANK FOUND
      QKR = QKR +1
      DPR(2) = DMX
      JP(5) = QMR
C- NOW TEST -IS THIS AN ITERATION TO DO ONLY 1 ROW AT A TIME/
      IF(QKR .EQ. QS) GO TO 480
      IF(JP(1) .EQ. IZ) GO TO 490
C( AT THIS POINT, EITHER STOP WITH ONE ROW OR TRY NEXT.
C HERE GET READY TO RETURN. RANK PARAMETER.
480   JP(3) = QKR
      RETURN
C IF ENOUGH TRIES TO D ALL ROWS PLUS 1 MORE, QUIT.
490   IF(QITR .EQ. QR1) GO TO 480
      QITR = QITR +1
      GO TO 200
  691 CALL LNCNT(1)
      WRITE (6,1691) QS,QNT
 1691 FORMAT  (' DIMENSION ERROR IN ANDRA NR=',I4,5X,'NR*NC='I4)
      RETURN
  692 CALL LNCNT(1)
      WRITE (6,1692)
 1692 FORMAT (' ERROR IN ANDRA, FINDS NO PIVOTS')
CHECK FOR DIAGONAL ALLOWING NO PIVOTS//
505   IF(JP(2) .EQ. IZ .OR. QKR .GT. QR1) GO TO 692
      GO TO 480
      END
ERASE SOURCE.DECOM$$
FTN DECOM$$,N,ISD=N
      SUBROUTINE DECOM (A,B,C,E,JL,DCM,KP,  D)
C- SUBR DECOM FINDS 3 MATRICES FROM WHICH USER CAN GET DECOMPOSITION
C INTO KRONECKER PRODUCT, POSSIBLY USING A SEPARATE PSEU CALL
C.. SUBR. DECOM  INFORMATION SYSTEMS CO. MAY 1969.  J W ANDREWS
      IMPLICIT REAL*8 (D), INTEGER*2 (Q)
      DOUBLE PRECISION A,B,C,E, D, PIV
      COMMON  /MAX/MAXRC
      DIMENSION A(400),B(400),C(400),E(400),  D(2000),
     2 JL(400),DCM(3),KP(4),NV(2),  ND(2)
      EQUIVALENCE (NV, I),  (PIV, ND(1)), (DFZO,IZ,QZ), (ND(2), J),
     1 (ND(1), QR1)
      DATA  DFZO /0.D0/
C- SET PARAMETERS TO CALL PSEU. TO FIND T TRANSFORMATION IN C.
C-- ASSUME INPUT A IS POS. DEFINITE SYMMETRIC
      QS = KP(3)
      QNT = QS*QS
C-ERR EXIT IF ABSURD DIMENSIONS.
      IF(QS .LT. 2 .OR. QNT .GT.MAXRC)  GO TO 691
C- SET COLUMN SIZE = RANK SIZE
      KP(4) = QS
      QL = KP(1)
C- SET SPECIAL PARAMS FOR PSEU CALL THESE ARE TO SUPPRESS THE WORK OF
C RE-INVERTING PSEUDO INVERSE IN THE CASE WHERE A SINGULAR...
      KP(1) =  - KP(1) -1
C- CALL PSEU P TO GET MATRIX T. IN C
CNOTE THE LAST 3 MATRICES OF THE 5 IN D USED ONLY IF PSEUP @ITERATES@
      I = KP(2)
          CALL PSEU P(A,B,C,E,DCM,KP,  D)
      KP(1) = QL
      IF(I .NE. IZ) GO TO  38
C/ PLEASE DO NOT TRY TO TAKE A.S.P. NAMES FOR MATRICES HERE.
C- SUCH MATRICES WERE NOT RETURNED BY ASP, NOR BY MY. ROUTINE.
C
      DO 13 I = 1, QNT
13    D(I) = C(I)
      NV(1)=QS
      NV(2)=QS
      ND(1) = 2
C- ND IS PART OF FLAG (PIV) RETURNED BY INV.
        CALL INV(D,NV,PIV,JL)
C- TEST TO SEE IF PIV IS ZERO = ERR, MATRIX NOT INVERTIBLE.  -ND1-
      IF(ND(1) .EQ. IZ) GO TO 692
38    CONTINUE
      KD = QNT
C- P IS TO HOLD PERMUTATION MATRIX SUCH THAT THAT  P*BE* PTR = ER
C- ER HAS ALL ONES MOVED TO EXTREME UPPER LEFT OF DIAGONAL.
C*NOW SET UP TO MAKE P  PERMUTATION MATRIX  P = D(KD +1)
      QR1 = QS +1
C- ZERO OUT P, WILL BE ZEROS AND ONES
      DO 39 I = 1, QNT
C-ZERO HOUSEKEEPING ARRAY  ONLY NEED FIRST COLUMN.
      JL(I) = IZ
      K = KD + I
39    D(K) = DFZO
      L =1
      M =1
      QL =1
      DO 780 K = 1, QS
      IF( B(L) ) 7803, 7801, 7803
7803  J = JL(QL)
CHECK FOR ROW OF DIAG. THAT NEEDS A 1 MOVED INTO IT
      IF( J .EQ. IZ) GO TO 786
      I = (K-1)*QS + J + KD
C-PUT 1 IN P TO MOVE K,K TO J,K POSITION  (2 PERMUTATIONS) P,, P TRANS
C*  THE EFFECT IS TO MOVE K,K TO J,K/ THENCE TO J,J.
      D(I) = 1.D0
      QL = QL +1
C/MARK THIS 1 AS ZRO TO BE FILLED-- IT IS MOVED UP AND OUTOF HERE
7801  JL(M) = K
      M = M +1
      GO TO 780
786   J = KD + L
C.MAKE PART OF IDENTITY AT 786 DON-T NEED TO MOVE 1 TO A HOLE.
      D(J) = 1.D0
780   L = QR1 + L
C RETURN. MATRICES COMPLETED E WITH IR 1@S DELIBERATELY LEFT OUT.
      RETURN
  691 CALL LNCNT(1)
      WRITE  (6,1691)  QS,QNT
 1691 FORMAT  (' DIMENSION ERROR IN DECOM  NC=',I4,5X,'NR*NC=',I4)
      RETURN
  692 CALL LNCNT(1)
      WRITE (6,1692)
 1692 FORMAT  (' ERROR IN DECOM   PIVOT=ZERO')
      KP(4)=-QS
      GO TO 38
      END
